--- title: "7-6 Shortness of breath" description: | Patient reported shortness of breath to day 28. author: - name: Rob Mahar affiliation: University of Melbourne - name: James Totterdell affiliation: University of Sydney date: today toc-depth: 4 --- # Preamble {#preamble} ```{r} #| label: pkgs #| code-summary: Load packages library (tidyverse)library (labelled)library (kableExtra)library (cmdstanr)library (posterior)library (bayestestR)library (bayesplot)library (matrixStats)library (patchwork)library (lubridate)library (ggdist)theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) + theme (panel.grid = element_blank (),strip.background = element_blank ()))bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) + theme (panel.grid = element_blank (),strip.background = element_blank ())))color_scheme_set ("red" )options (digits = 4 )``` ```{r} #| label: data #| code-summary: Prepare dataset :: load_all ()<- read_all_no_daily ()# Anticoagulation set (intention-to-treat, including missing) <- all_dat %>% filter_acs_itt () %>% transmute_model_cols_grp_aus_nz ()# Anticoagulation set (intention-to-treat, excluding missing) <- acs_itt_dat %>% filter (! is.na (out_sob)) # Concurrent enrolments for C2 <- acs_itt_dat %>% filter_concurrent_intermediate ()<- acs_itt_concurc2_dat %>% filter (! is.na (out_sob))# Concurrent enrolments for C3 <- acs_itt_dat %>% filter_concurrent_std_aspirin ()<- acs_itt_concurc3_dat %>% filter (! is.na (out_sob))# Concurrent enrolments for C4 <- acs_itt_dat %>% filter_concurrent_therapeutic ()<- acs_itt_concurc4_dat %>% filter (! is.na (out_sob))``` ```{r} #| label: models #| code-summary: Load the Stan models <- cmdstan_model (file.path ("stan" , "binary" , "logistic_site_epoch.stan" ))<- cmdstan_model (file.path ("stan" , "binary" , "logistic_site.stan" ))<- cmdstan_model (file.path ("stan" , "binary" , "logistic_epoch.stan" )) <- cmdstan_model (file.path ("stan" , "binary" , "logistic.stan" ))``` ```{r} #| label: helpers #| code-summary: Helper functions <- function (dat, format = "html" ) {<- dat %>% mutate (CAssignment = factor (levels = c ("C1" , "C2" , "C3" , "C4" ),labels = str_replace (intervention_labels ()$ CAssignment[- 1 ], "<br>" , " " ))) %>% group_by (CAssignment) %>% summarise (Patients = n (),Known = sprintf ("%i (%.1f)" , sum (! is.na (out_sob)), 100 * mean (! is.na (out_sob))),Missing = sprintf ("%i (%.1f)" , sum (is.na (out_sob)), 100 * mean (is.na (out_sob))),` Shortness of breath day 28 ` = sprintf ("%i (%.1f)" , sum (out_sob, na.rm = TRUE ), 100 * mean (out_sob, na.rm = TRUE )),%>% bind_rows (%>% group_by (CAssignment = "Overall" ) %>% summarise (Patients = n (),Known = sprintf ("%i (%.1f)" , sum (! is.na (out_sob)), 100 * mean (! is.na (out_sob))),Missing = sprintf ("%i (%.1f)" , sum (is.na (out_sob)), 100 * mean (is.na (out_sob))),` Shortness of breath day 28 ` = sprintf ("%i (%.1f)" , sum (out_sob, na.rm = TRUE ), 100 * mean (out_sob, na.rm = TRUE )),%>% mutate (CAssignment = fct_inorder (CAssignment)) %>% gather (key, value, - CAssignment, factor_key = T) %>% spread (key, value)colnames (tab)[1 ] <- "n ( \\ %)" if (format == "latex" ) {colnames (tab) <- linebreak (colnames (tab), align = "c" , linebreaker = "<br>" )kable (tab,format = format,align = "lrrrrr" ,escape = FALSE ,booktabs = TRUE %>% kable_styling (font_size = 10 , latex_options = "HOLD_position" ) %>% row_spec (nrow (tab), bold = TRUE )<- function (vars = NULL , outcome = NULL , beta_sd = NULL , ctr = contr.orthonorm) {<- make_X_design (dat, vars, ctr)<- sum (grepl ("rand" , colnames (X))) - 1 <- c (2.5 , rep (1 , nXassign), beta_sd)<- dat[[outcome]]<- dim (X)[1 ]<- dim (X)[2 ] <- dat$ epoch<- max (dat$ epoch)<- dat[["ctry_num" ]]<- max (region)<- dat[["site_num" ]]<- max (site)<- region_by_site <- dat %>% :: count (ctry_num, site_num) %>% pull (ctry_num)list (N = N, K = K, X = X, y = y,M_region = M_region, region = region,M_site = M_site, site = site,M_epoch = M_epoch, epoch = epoch,region_by_site = region_by_site,beta_sd = beta_sd)<- function (vars = NULL , outcome = NULL , beta_sd = NULL , ctr = contr.orthonorm) {<- make_X_design (dat, vars, ctr, includeA = FALSE )<- sum (grepl ("rand" , colnames (X))) - 1 <- c (2.5 , rep (1 , nXassign), beta_sd)<- dat[[outcome]]<- dim (X)[1 ]<- dim (X)[2 ] list (N = N, K = K, X = X, y = y, beta_sd = beta_sd)``` # Descriptive ## Anticoagulation ```{r} #| label: tbl-anticoagulation-sob #| code-summary: Anticoagulation #| tbl-cap: Summary of outcome by anticoagulation intervention. save_tex_table (make_summary_table (acs_itt_dat, "latex" ),"outcomes/secondary/7-6-anticoagulation-summary" )make_summary_table (acs_itt_dat)``` ## Age ```{r} #| label: fig-age-sob #| code-summary: Relationship age to shortness of breath #| fig-cap: | #| Relationship (logistic regression linear in age) #| between age at entry and shortness of breath at day 28. <- acs_itt_dat %>% :: count (out_sob, AgeAtEntry) %>% spread (out_sob, n, fill = 0 ) %>% mutate (n = ` 0 ` + ` 1 ` + ` <NA> ` ,p = ` 1 ` / (` 1 ` + ` 0 ` ))<- glm (cbind (` 1 ` , ` 0 ` ) ~ AgeAtEntry, data = agedat, family = binomial ())<- agedat %>% mutate (ypred = predict (agemod, newdata = agedat, type = "response" )<- ggplot (agedat, aes (AgeAtEntry, n)) + geom_col (colour = "grey40" , fill = "grey40" ) + geom_vline (xintercept = 60 , linetype = 2 ) + labs (y = "Number of \n participants" ,x = "Age at entry" )<- ggplot (agedat, aes (AgeAtEntry, p)) + geom_point () + geom_vline (xintercept = 60 , linetype = 2 ) + geom_line (aes (y = ypred)) + labs (y = "Proportion \n shortness of breath \n day 28" , x = "Age at entry" )<- p1 | p2<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-6-age.pdf" ), p, height = 2 , width = 6 )``` ## Country ```{r} #| label: fig-country-sob #| code-summary: Relationship country to outcome #| fig-cap: Shortness of breath by country. <- acs_itt_dat %>% :: count (Country = factor (levels = c ("IN" , "AU" , "NP" , "NZ" ),labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )), out_sob) %>% group_by (Country) %>% spread (out_sob, n, fill = 0 ) %>% mutate (n = ` 1 ` + ` 0 ` + ` <NA> ` ,p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )<- ggplot (tdat, aes (Country, n)) + geom_col () + labs (y = "Number of \n participants" , x = "Country of enrolment" )<- ggplot (tdat, aes (Country, p_1)) + geom_point () + labs (y = "Proportion \n shortness of breath \n day 28" , x = "Country of enrolment" ) + ylim (0 , NA )<- p1 | p2<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-6-country.pdf" ), p, height = 2 , width = 6 )``` ## Site ```{r} #| label: fig-site-7-6 #| code-summary: Relationship site to outcome #| fig-cap: Shortness of breath at day 28 by site within country. <- all_dat %>% filter_acs_itt () %>% :: count (Country_lab = Country,Site_lab = fct_infreq (Location),Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" )),Site = PT_LocationName,%>% group_by (Country, Site) %>% spread (out_sob, n, fill = 0 ) %>% mutate (n = ` 1 ` + ` 0 ` + ` <NA> ` ,p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )%>% ungroup ()<- ggplot (tdat, aes (Site_lab, n)) + facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) + geom_col () + labs (y = "Number of \n participants" , x = "" ) + theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),panel.border = element_rect (fill = NA ))<- ggplot (tdat, aes (Site_lab, p_1)) + facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) + geom_point () + labs (y = "Proportion \n shortness of breath \n day 28" , x = "Site of enrolment" ) + theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),panel.border = element_rect (fill = NA ))<- p1 / p2<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-6-country-site.pdf" ), p, height = 4 , width = 6.25 )``` ## Calendar Time ```{r} #| label: fig-cal-po #| code-summary: Relationship calendar date to outcome #| fig-cap: | #| Relationship between calendar date and the primary outcome. <- all_dat %>% filter_acs_itt () %>% :: count (out_sob, yr = year (RandDate), mth = month (RandDate)) %>% spread (out_sob, n, fill = 0 ) %>% mutate (p = ` 1 ` / (` 1 ` + ` 0 ` ),n = ` 1 ` + ` 0 ` + ` <NA> ` )<- ggplot (caldat, aes (mth, p)) + facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) + geom_point () + labs (y = "Proportion \n shortness of breath \n day 28" , x = "Calendar date (month of year)" ) + scale_x_continuous (breaks = 1 : 12 ) + ylim (0 , NA )<- ggplot (caldat, aes (mth, n)) + facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) + geom_col () + labs (y = "Number of \n participants" , x = "Calendar date (month of year)" ) + scale_x_continuous (breaks = 1 : 12 )<- p2 | p1<- file.path ("outputs" , "figures" , "outcomes" , "secondary" )ggsave (file.path (path, "7-6-calendar-time.pdf" ), p, height = 2 , width = 6 )``` # Modelling ## Primary model ```{r} #| label: model1-sampling #| code-summary: Data and sampling <- 36127 <- make_stan_data (dat = acs_itt_nona_dat, vars = c ("inelgc3" , "agegte60" , "ctry" ),outcome = "out_sob" ,beta_sd = c (10 , 2.5 , 1 , 1 ))<- capture.output (<- model_full$ sample (data = mdat,seed = seed,refresh = 0 ,iter_warmup = 1000 ,iter_sampling = 2500 ,chains = 8 ,parallel_chains = 8 ,adapt_delta = 0.95 $ save_object (file.path ("outputs" , "models" , "secondary" , "7-6.rds" ))<- as_draws_rvars (mfit$ draws ())names (mdrws$ beta) <- colnames (mdat$ X)names (mdrws$ gamma_epoch) <- acs_itt_nona_dat %>% dplyr:: count (epoch, epoch_lab) %>% pull (epoch_lab)names (mdrws$ gamma_site) <- acs_itt_nona_dat %>% dplyr:: count (site_num, site) %>% pull (site)<- factor (mdat$ region_by_site, labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))$ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]$ Atrt <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),"Therapeutic vs low" = exp (mdrws$ Ctrt[3 ]),"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]),"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[3 ]),"Low with aspirin vs therapeutic" = exp (mdrws$ Ctrt[2 ] - mdrws$ Ctrt[3 ])$ OR <- c (setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),setNames (exp (mdrws$ beta[7 : 8 ]), c ("Ineligible aspirin" , "Age \u2265 60" )),setNames (exp (mdrws$ beta[9 : 10 ]), c ("Australia/New Zealand" , "Nepal" ))``` ```{r} #| label: odds-ratio-summary-table-primary-model-save #| code-summary: Odds ratio summary table #| tbl-cap: Posterior summaries for model parameters (fixed-effects). save_tex_table (odds_ratio_summary_table (mdrws$ OR, "latex" ),"outcomes/secondary/7-6-primary-model-acs-itt-summary-table" )odds_ratio_summary_table (mdrws$ OR)``` ```{r} #| code-summary: Odds ratio summary for epoch and site #| fig-cap: Summary of epoch and site posterior odds ratios. <- plot_epoch_site_terms (exp (mdrws$ gamma_epoch),exp (mdrws$ gamma_site),<- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-6-acs-itt-epoch-site-terms.pdf" )ggsave (pth, p, width = 6 , height = 4.5 )``` ```{r} #| fig-cap: Posterior densities for treatment comparisons. <- plot_or_densities (mdrws$ compare)``` ### Diagnostics ```{r} $ summary (variables = c ("beta" , "gamma_site" , "gamma_epoch" ))$ diagnostic_summary ()``` ### Trace plots ```{r} mcmc_trace (mdrws["beta" ])mcmc_trace (mdrws["gamma_site" ])mcmc_trace (mdrws["gamma_epoch" ])mcmc_trace (mdrws[c ("tau_site" , "tau_epoch" )])``` ### Posterior predictive ```{r} #| label: primary-acs-itt-ppc #| fig-cap: Posterior predictive distribution versus observed proportion (red diamond). <- mdrws$ y_ppc<- bind_cols (acs_itt_nona_dat, tibble (y_ppc = y_ppc))<- function (grp) {%>% group_by (grp = !! grp) %>% summarise (mean_y = mean (out_sob),rvar_mean_y_ppc = rvar_mean (y_ppc)<- function (dat, lab = "" ) {%>% ggplot (aes (y = grp, xdist = rvar_mean_y_ppc)) + stat_interval (size = 2 ) + geom_point (aes (x = mean_y, y = 1 : nrow (dat)), colour = "red" , shape = 23 ) + labs (x = "Posterior predictive \n proportion" , y = lab) <- grp_ppc (sym ("CAssignment" ))<- grp_ppc (sym ("country" ))<- grp_ppc (sym ("epoch" ))<- ppc_dat %>% group_by (Country = country, grp = site_raw) %>% summarise (mean_y = mean (out_sob),rvar_mean_y_ppc = rvar_mean (y_ppc)<- plot_grp_ppc (ppc_C, "Anticoagulation \n intervention" ) + labs (x = "" )<- plot_grp_ppc (ppc_ctry, "Country" ) + labs (x = "" )<- plot_grp_ppc (ppc_epoch, "Epoch" ) + labs (x = "" )<- plot_grp_ppc (ppc_site %>% filter (Country == "IN" ), "Sites \n India" ) + labs (x = "" )<- plot_grp_ppc (ppc_site %>% filter (Country == "AU" ), "Sites \n Australia" ) + labs (x = "" )<- plot_grp_ppc (ppc_site %>% filter (Country == "NP" ), "Sites \n Nepal" )<- plot_grp_ppc (ppc_site %>% filter (Country == "NZ" ), "Sites \n NZ" )<- ((p3 | p1 / p2) / | p5) / (p6 | p7) + plot_layout (heights = c (3 , 1 )) ) ) + plot_layout (heights = c (1 , 1.5 ), guides = "collect" )<- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,"7-6-primary-model-acs-itt-ppc.pdf" )ggsave (pth, p, width = 6 , height = 5.5 )``` ## Sensitivity: Concurrent Controls ### Intermediate ```{r} #| code-summary: Overview of outcome #| tbl-cap: Summary of outcome for intermediate dose concurrent enrolments. save_tex_table (make_summary_table (acs_itt_concurc2_dat, "latex" ),"outcomes/secondary/7-6-anticoagulation-concurrent-intermediate-summary" )make_summary_table (acs_itt_concurc2_dat)``` ```{r} #| code-summary: Fit model <- make_stan_data_concurrent (c ("agegte60" , "ctry" ),"out_sob" ,c (2.5 , 1 , 1 )<- capture.output (<- model_fixed_only$ sample (data = mdat,seed = 38135 ,chains = 8 ,parallel_chains = 8 ,iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 <- as_draws_rvars (mfit$ draws (c ("beta" )))names (mdrws$ beta) <- colnames (mdat$ X)$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 ]$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]))$ OR <- c (setNames (exp (mdrws$ Ctrt), c ("Intermediate" )),setNames (exp (mdrws$ beta[- (1 : 2 )]), c ("Age \u2265 60" , "Australia/New Zealand" , "Nepal" ))``` ```{r} #| code-summary: Odds ratio summary #| tbl-cap: Odds ratio summary table for model. save_tex_table (odds_ratio_summary_table (mdrws$ OR, "latex" ),"outcomes/secondary/7-6-primary-model-acs-itt-concurrent-intermediate-summary-table" )odds_ratio_summary_table (mdrws$ OR)``` ```{r} #| code-summary: Posterior contrast #| tbl-cap: Posterior densities for treatment comparisons, ACS-ITT-intermediate <- plot_or_densities (mdrws$ compare)``` #### Diagnostics ```{r} $ summary ()$ diagnostic_summary ()``` #### Trace plots ```{r} mcmc_trace (mdrws["beta" ])``` ### Low-dose with aspirin ```{r} #| code-summary: Overview of outcome #| tbl-cap: Summary of outcome for low-dose with aspirin concurrent enrolments. save_tex_table (make_summary_table (acs_itt_concurc3_dat, "latex" ),"outcomes/secondary/7-6-anticoagulation-concurrent-stdaspirin-summary" )make_summary_table (acs_itt_concurc3_dat)``` ```{r} #| code-summary: Fit model <- make_stan_data_concurrent (c ("agegte60" , "ctry" ),"out_sob" ,c (2.5 , 1 )<- capture.output (<- model_fixed_only$ sample (data = mdat,seed = 13578 ,chains = 8 ,parallel_chains = 8 ,iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 <- as_draws_rvars (mfit$ draws (c ("beta" )))names (mdrws$ beta) <- colnames (mdat$ X)$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 3 ]$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))$ OR <- c (setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" )),setNames (exp (mdrws$ beta[- (1 : 3 )]), c ("Age \u2265 60" , "Australia/New Zealand" ))``` ```{r} #| code-summary: Odds ratio summary #| tbl-cap: Odds ratio summary table for model. save_tex_table (odds_ratio_summary_table (mdrws$ OR, "latex" ),"outcomes/secondary/7-6-primary-model-acs-itt-concurrent-stdaspirin-summary-table" )odds_ratio_summary_table (mdrws$ OR)``` ```{r} #| code-summary: Posterior contrast #| tbl-cap: Posterior densities for treatment comparisons, ACS-ITT-aspirin <- plot_or_densities (mdrws$ compare)``` #### Diagnostics ```{r} $ summary ()$ diagnostic_summary ()``` #### Trace plots ```{r} mcmc_trace (mdrws["beta" ])``` ### Therapeutic ```{r} #| code-summary: Overview of outcome #| tbl-cap: Summary of outcome for therapeutic dose concurrent enrolments. save_tex_table (make_summary_table (acs_itt_concurc4_dat, "latex" ),"outcomes/secondary/7-6-anticoagulation-concurrent-therapeutic-summary" )make_summary_table (acs_itt_concurc4_dat)``` ```{r} #| code-summary: Fit model <- make_stan_data_concurrent (c ("agegte60" , "ctry" ),"out_sob" ,c (2.5 , 1 , 1 )<- capture.output (<- model_fixed_only$ sample (data = mdat,seed = 13578 ,chains = 8 ,parallel_chains = 8 ,iter_warmup = 1000 ,iter_sampling = 2500 ,refresh = 0 <- as_draws_rvars (mfit$ draws (c ("beta" )))names (mdrws$ beta) <- colnames (mdat$ X)$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 3 ]$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),"Therapeutic vs low" = exp (mdrws$ Ctrt[2 ]),"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))$ OR <- c (setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Therapeutic" )),setNames (exp (mdrws$ beta[- (1 : 3 )]), c ("Age \u2265 60" , "India" , "Australia/New Zealand" ))``` ```{r} #| code-summary: Odds ratio summary #| tbl-cap: Odds ratio summary table for model. save_tex_table (odds_ratio_summary_table (mdrws$ OR, "latex" ),"outcomes/secondary/7-6-primary-model-acs-itt-concurrent-therapeutic-summary-table" )odds_ratio_summary_table (mdrws$ OR)``` ```{r} #| code-summary: Posterior contrast #| tbl-cap: Posterior densities for treatment comparisons, ACS-ITT-therapeutic <- plot_or_densities (mdrws$ compare)``` #### Diagnostics ```{r} $ summary ()$ diagnostic_summary ()``` #### Trace plots ```{r} mcmc_trace (mdrws["beta" ])``` 
